#####
#Replication Code for Jesse Driscoll and F. Daniel Hidalgo (2014): "Intended and Unintended Consequences of Democracy Promotion Assistance to Georgia After The Rose Revolution"
#Research and Politics
#####


###Set working directory to directory containing data
setwd('/Users/dhidalgo/Dropbox/projects/georgia_experiment/rcode/replication_code/')

###
load('hidalgo_driscoll_2014_replication.RData')
###

##Following packages are required
library(apsrtable)
library(sandwich)
library(survey)
library(ggplot2)
library(gdata)
library(stargazer)
library(plyr)
library(Hmisc)

#####
#Table 1
#####

#adjust presentation details of apsrtable
invisible(setMethod("modelInfo", "summary.lm", function(x) {
  env <- parent.frame()
  digits <- evalq(digits,env)
  model.info <- list(
    "Covariates" = (
      ifelse(!is.na(charmatch("voters.4.17",rownames(coef(x)))),
             "\\textnormal{  Yes}", "\\textnormal{  No}")),
    "N"= (formatC(sum(x$df[1:2]),format="d")))
  class(model.info) <- "model.info"
  return(model.info)
} ))


covariates <- c("block", "voters.4.17", "turnout2006", "pctPresParty2006")
lm.complaints.dummy <- (lm(complaints.dummy~tr.complaints + block, data=exp.data))
lm.complaints.dummy$se <- vcovHC(lm.complaints.dummy,type = "HC0")
lm.complaints.dummy.cov <- (lm(paste("complaints.dummy~tr.complaints+block",paste(covariates, collapse="+"), sep = "+"), data=exp.data))
lm.complaints.dummy.cov$se <- vcovHC(lm.complaints.dummy.cov,type = "HC0")
lm.total.complaints <- (lm(total.complaints~tr.complaints + block, data=exp.data))
lm.total.complaints$se <- vcovHC(lm.total.complaints,type = "HC0")
lm.total.complaints.cov <- (lm(paste("total.complaints~tr.complaints + block",paste(covariates, collapse="+"), sep = "+"), data=exp.data))
lm.total.complaints.cov$se <- vcovHC(lm.total.complaints.cov,type = "HC0")


complaints.caption <- "The average treatment effect of the complaints treatment on whether or not any complaint was made and the  number of complaints in the 2008 elections. For covariate adjusted results, the adjustment variables are number of registered voters, 2006 turnout, and 2006 vote share of the ruling party. Standard errors are heteroskedasticity consistent and account for blocking."

complaints_table <- apsrtable(lm.complaints.dummy,lm.complaints.dummy.cov, lm.total.complaints, lm.total.complaints.cov, Sweave = FALSE, omitcoef = expression(grep(pattern=paste("Intercept", paste(covariates,collapse="|"), sep = "|"),x=coefnames)), coef.names = c("Complaints Treatment"), stars = "default", model.names = c("1","2","3","4"), caption = complaints.caption, caption.position = "above", label = "tab:dv.complaints")
complaints_table[2] <- paste("\n \n  & \\multicolumn{ 2 }{ c }{ Complaints Dummy } & \\multicolumn{ 2 }{ c }{ Total Complaints }\\\\ \n", complaints_table[2], sep=" ")
#complaints_table[1] <- gsub(x = complaints_table[1], pattern = "\n \n", replacement = "\n \\\\centering \n")
complaints_table[1] <- gsub(x = complaints_table[1], pattern = "\\!ht\\]\n", replacement = "\\!ht\\]\n \\\\centering \n")
writeLines(complaints_table, con = "./complaints_table.tex")


####
#Table 2
###


ols.turnout <- (lm((turnout08)~tr.complaints*tr.list+block, data=exp.data))
ols.turnout$se <- vcovHC(ols.turnout,type = "HC0")
ols.turnout.cov <- (lm(paste("turnout08~tr.complaints*tr.list+block",paste(covariates, collapse="+"), sep = "+"), data=exp.data))
ols.turnout.cov$se <- vcovHC(ols.turnout.cov,type = "HC0")

ols.pres <- (lm((pctPresParty08)~tr.complaints*tr.list+block,data=exp.data))
ols.pres$se <- vcovHC(ols.pres,type = "HC0")
ols.pres.cov <- (lm(paste("pctPresParty08~tr.complaints*tr.list+block",paste(covariates, collapse="+"), sep = "+"),data=exp.data))
ols.pres.cov$se <- vcovHC(ols.pres.cov,type = "HC0")

unanticipated.caption <- "The average treatment effect of the complaints treatment on turnout and the vote share of the ruling party. For covariate adjusted results, the adjustment variables   are number of registered voters, 2006 turnout, and 2006 vote share   of the ruling party. Standard errors are heteroskedasticity consistent and account for blocking."

unanticipated_table <- apsrtable(ols.turnout, ols.turnout.cov, ols.pres, ols.pres.cov, Sweave = FALSE, omitcoef = expression(grep(pattern=paste("Intercept", "tr.list", paste(covariates,collapse="|"), sep = "|"),x=coefnames)), coef.names = c("Complaints Treatment"), model.names = c("1","2","3","4"), caption = unanticipated.caption, caption.position = "above", label = "tab:dv.unanticipated")
unanticipated_table[2] <- paste("\n \n  & \\multicolumn{ 2 }{ c }{ Turnout (\\%) } & \\multicolumn{ 2 }{ c }{ Pres. Vote (\\%) }\\\\ \n", unanticipated_table[2], sep=" ")
#unanticipated_table[1] <- gsub(x = unanticipated_table[1], pattern = "\n \n", replacement = "\n \\\\centering \n")
unanticipated_table[1] <- gsub(x = unanticipated_table[1], pattern = "\\!ht\\]\n", replacement = "\\!ht\\]\n \\\\centering \n")
writeLines(unanticipated_table, con = "unanticipated_table.tex")



### Table 3

exp.design.08 <- svydesign(ids = ~psu, weight = ~wt, data = svy.asexp.data.08) 
svy.asexp.turnout08 <- (svyglm(turnout08 ~ treat + regvoters.mos, design = exp.design.08))
svy.asexp.pres08 <- (svyglm(pres.pct.08 ~treat + regvoters.mos, design = exp.design.08))

exp.design.10 <- svydesign(ids = ~psu.id, weights = ~wt, strata = ~strata, data = svy.asexp.data.10)
svy.asexp.majturn10 <- (svyglm(maj.turnout.pct~treat + numvoters.frame, design = exp.design.10))
svy.asexp.pres10 <- (svyglm(presparty10.pct~treat  + numvoters.frame, design = exp.design.10))

conf.int <- function(x){
  paste("(", round(x$coef[2] - 1.96  * SE(x)[2], 2), ", ",
        round(x$coef[2] + 1.96 * SE(x)[2], 2),")", sep = "")     
}

svy.asexp.table <- data.frame(matrix(ncol = 4, nrow = 0))
svy.asexp.table <- rbind(svy.asexp.table, round(c(coef(svy.asexp.turnout08)[2], coef(svy.asexp.majturn10)[2], coef(svy.asexp.pres08)[2], coef(svy.asexp.pres10)[2]), 2),
                         round(c(SE(svy.asexp.turnout08)[2], SE(svy.asexp.majturn10)[2], SE(svy.asexp.pres08)[2], SE(svy.asexp.pres10)[2]), 2),
                         c(conf.int(svy.asexp.turnout08), conf.int(svy.asexp.majturn10), conf.int(svy.asexp.pres08), conf.int(svy.asexp.pres10)),
                         as.character(c(nrow(svy.asexp.turnout08$data), nrow(svy.asexp.majturn10$data), nrow(svy.asexp.pres08$data), nrow(svy.asexp.pres10$data)))
)
names(svy.asexp.table) <- c("2008", "2010", "2008", "2010")

svy.asexp.caption <- c("This table shows the estimated effects of a precinct being included in two different political surveys. The first survey was conducted immediately before the 2008 parliamentary election. The second survey was carried out before the 2010 local elections. Treatment effect estimates account for unequal probability of selection and control for precinct size, the variable used to sample precincts into the ``treatment''.")

latex(svy.asexp.table, cdec = c(2,2,2,2), dcolumn = FALSE, ctable = TRUE,
      cgroup = c("Turnout \\%", "Pres. Party Vote \\%"), n.cgroup = c(2,2),
      rowname = c("Estimate", "SE", "95\\% CI", "N"), rowlabel = "",
      col.just = rep("c", 4), insert.bottom = svy.asexp.caption, where = "htp",
      label = "svy.asexp.table", caption = "Survey ``as an'' Experiment",
      file = "svy.asexp.table.tex")


###
## Figure 1
##

hetero.job.pressure.hi.mod <- (lm(turnout08~tr.complaints.only  + turnout2006+pctPresParty2006+num.reg.voters + block, data = exp.data, subset = inter.job.pressure.dummy==1))
hetero.job.pressure.lo.mod <- (lm(turnout08~tr.complaints.only + turnout2006+pctPresParty2006+num.reg.voters + block, data = exp.data, subset = inter.job.pressure.dummy==0))
hetero.urban.hi.mod <- (lm(turnout08~tr.complaints.only + turnout2006+pctPresParty2006+num.reg.voters  + block, data = exp.data, subset = urban==1))
hetero.urban.lo.mod <- (lm(turnout08~tr.complaints.only +turnout2006+pctPresParty2006+num.reg.voters  + block, data = exp.data, subset = urban==0))
hetero.empl.camp.hi.mod <- (lm(turnout08~tr.complaints.only + turnout2006+pctPresParty2006+num.reg.voters + block, data = exp.data, subset = inter.empl.camp.dummy==1))
hetero.empl.camp.lo.mod <- (lm(turnout08~tr.complaints.only + turnout2006+pctPresParty2006+num.reg.voters  + block, data = exp.data, subset = inter.empl.camp.dummy==0))

hetero.exp.table <- data.frame(Estimate = c(coef(hetero.job.pressure.hi.mod)[2],
                                            coef(hetero.job.pressure.lo.mod)[2],
                                            coef(hetero.empl.camp.hi.mod)[2],
                                            coef(hetero.empl.camp.lo.mod)[2],
                                            coef(hetero.urban.hi.mod)[2],
                                            coef(hetero.urban.lo.mod)[2]),
                               SE = c(SE(hetero.job.pressure.hi.mod)[2],
                                      SE(hetero.job.pressure.lo.mod)[2],
                                      SE(hetero.empl.camp.hi.mod)[2],
                                      SE(hetero.empl.camp.lo.mod)[2],
                                      SE(hetero.urban.hi.mod)[2],
                                      SE(hetero.urban.lo.mod)[2])
)
hetero.exp.table$Variable <- c('Employer Pressure to Vote','Employer Pressure to Vote', 'Employer Pressure to Campaign', 'Employer Pressure to Campaign', 'Rural or Urban', 'Rural or Urban')
hetero.exp.table$Strata <- c('Above Median', 'Below Median', 'Above Median', 'Below Median', 'Urban', 'Rural')

ggplot(hetero.exp.table, aes(x = Strata, y = Estimate)) + geom_pointrange(size = 1.5, aes(ymax = Estimate + 1.96 * SE, ymin = Estimate - 1.96 * SE)) + facet_wrap( ~ Variable, nrow = 1, scales = "free_x") + theme_bw() + geom_hline(yintercept = 0, lty = 2) + ylab("Treatment Effect on Turnout (%)")


##
###Table 4
##

block.ate <- function(dv, treat, block){
  data <- data.frame(dv,treat,block)
  n <-  nrow(data)
  diff.means <- ddply(data, "block", summarise, treat.eff = mean(dv[treat==1]) - mean(dv[treat==0]), se = sqrt((var(dv[treat==1])/sum(treat)) + (var(dv[treat==0])/sum(1-treat))),  wt = length(treat)/n)
  ate <- weighted.mean(diff.means$treat.eff, diff.means$wt)
  se.ate <-sqrt(with(diff.means, mean(se^2 * wt^2/sum(wt^2))))
  list(ate = ate, se = se.ate, n = length(treat))
}

ate.refuse.dk <- block.ate(dv = svy10.analysis.data$sum.refuse.dk, treat = svy10.analysis.data$treat, block = svy10.analysis.data$psu)
ate.refuse <- block.ate(dv = svy10.analysis.data$sum.refuse, treat = svy10.analysis.data$treat, block = svy10.analysis.data$psu)
ate.dk <- block.ate(dv = svy10.analysis.data$sum.dontknow, treat = svy10.analysis.data$treat, block = svy10.analysis.data$psu)


####
#Table 5
###


ate.democracy.low <- with(svy10.analysis.data[svy10.analysis.data$psu.is.georgia.democracy<median(svy10.analysis.data$psu.is.georgia.democracy),], block.ate(dv = sum.refuse.dk, treat = treatment, block = psu))
ate.democracy.hi <- with(svy10.analysis.data[svy10.analysis.data$psu.is.georgia.democracy>median(svy10.analysis.data$psu.is.georgia.democracy),], block.ate(dv = sum.refuse.dk, treat = treatment, block = psu))

ate.clean.low <- with(svy10.analysis.data[svy10.analysis.data$psu.local.election.clean<median(svy10.analysis.data$psu.local.election.clean),], block.ate(dv = sum.refuse.dk, treat = treatment, block = psu))
ate.clean.hi <- with(svy10.analysis.data[svy10.analysis.data$psu.local.election.clean>median(svy10.analysis.data$psu.local.election.clean),], block.ate(dv = sum.refuse.dk, treat = treatment, block = psu))

ate.unm.low <- with(svy10.analysis.data[svy10.analysis.data$psu.UNM.closest<median(svy10.analysis.data$psu.UNM.closest),], block.ate(dv = sum.refuse.dk, treat = treatment, block = psu))
ate.unm.hi <- with(svy10.analysis.data[svy10.analysis.data$psu.UNM.closest>median(svy10.analysis.data$psu.UNM.closest),], block.ate(dv = sum.refuse.dk, treat = treatment, block = psu))


ate.democracy.low
ate.democracy.hi
ate.clean.low
ate.clean.hi
ate.unm.low
ate.unm.hi

### Table 6
bal.var <- c('voters.4.17', 'turnout2006', 'pctPresParty2006', 'pre.nonresp', 'pre.job.pressure', 'pre.not.democr', 'pre.ballot.stuff', 'off.bias', 'harass.voters', 'harass.act')

expBalStat <- function(dv, data = exp.data){
  data$dv <- data[, dv]
  complaints.model <- (lm(dv ~ tr.complaints + block, data = data))
  voterlist.model <- lm(dv ~ tr.list + block, data = data, subset = urban == 1)
  list(complaints_coef = coef(complaints.model)['tr.complaints'], se_coef = sqrt(vcovHC(complaints.model,type = "HC0")[2,2]),
       voterlist_coef = coef(voterlist.model)['tr.list'], se_coef = sqrt(vcovHC(voterlist.model,type = "HC0")[2,2]))
}


bal_stats <- data.frame(sapply(bal.var, expBalStat))
variables <- c('Registered Voters (2007)', 'Turnout (2006)', 'Presidential Party Vote Pct. (2006)', 'Pre-Treatment Non-Response Pct.', 'Pct. Saying Employer Pressure is Widespread', 'Pct. Believing Georgia Not a Democracy', 'Pct. Believing Ballot Stuffing is Widespread', 'Pct. Believing Officals Biased to One Party', 'Pct. Believing Voter Harassment is Widespread', 'Pct. Believing Activist Harassment is Widespread')
complaints.bal.coefs <- round(unlist(bal_stats[1,]), 2)
complaints.bal.se <- sapply(round(unlist(bal_stats[2,]), 2), function(x) paste("(", x, ")", sep =""))
voterlist.bal.coefs <- round(unlist(bal_stats[3,]), 2)
voterlist.bal.se <- sapply(round(unlist(bal_stats[4,]), 2), function(x) paste("(", x, ")", sep =""))

bal_stats_table <- data.frame(
  Covariate = matrix(interleave(variables, rep("", length(variables))),ncol = 1), 
  'Complaints Treatment' = matrix(interleave(complaints.bal.coefs, complaints.bal.se),ncol = 1 ),
  'Voter List Treatment' = matrix(interleave(voterlist.bal.coefs, voterlist.bal.se),ncol =1)
)
bal_stats_caption <- "Table shows covariate balance on selected pre-treatment covariates for both the complaints intervention (N=84) and the voter list intervention (n = 48). Point estimates are from a regression with a treatment dummy and  block fixed effects. Standard errors (in parentheses) are heteroskedasticity consistent."

latex(bal_stats_table, cdec = c(2,2,2), dcolumn = FALSE, ctable = TRUE,
      rowlabel = "", rowname = NULL, colheads = c('Covariate', "Complaints Treatment", "Voter List Treatment"),
      col.just = c("l", "c", "c"), insert.bottom = bal_stats_caption, where = "htp",
      label = "bal.stats.table", caption = "Balance Statistics for the 2008 Field Experiment",
      file = "exp_bal_stats.tex")

##Table 7

svy.asexp.past_turnout08 <- (svyglm(turnout.06 ~ treat, design = exp.design.08))
svy.asexp.past_pres08 <- (svyglm(pct.presparty.06 ~ treat, design = exp.design.08))
svy.asexp.past_regvoters08 <- (svyglm(reg.voters.06 ~ treat, design = exp.design.08))

svy.asexp.past_turnout10 <- (svyglm(turnout08 ~treat + numvoters.frame, design = exp.design.10))
svy.asexp.past_pres10 <- (svyglm(pctPresParty08 ~treat + numvoters.frame, design = exp.design.10))
svy.asexp.past_regvoters10 <- (svyglm(regvoters.08 ~treat + numvoters.frame, design = exp.design.10))



svy.asexp.balance.table <- data.frame(matrix(ncol = 4, nrow = 0))
svy.asexp.balance.table <- rbind(svy.asexp.balance.table, round(c(coef(svy.asexp.past_turnout08)[2], coef(svy.asexp.past_turnout10)[2], coef(svy.asexp.past_pres08)[2], coef(svy.asexp.past_pres10)[2]), 2),
                                 round(c(SE(svy.asexp.past_turnout08)[2], SE(svy.asexp.past_turnout10)[2], SE(svy.asexp.past_pres08)[2], SE(svy.asexp.past_pres10)[2]), 2),
                                 c(conf.int(svy.asexp.past_turnout08), conf.int(svy.asexp.past_turnout10), conf.int(svy.asexp.past_pres08), conf.int(svy.asexp.past_pres10)),
                                 as.character(c(nrow(svy.asexp.turnout08$data), nrow(svy.asexp.majturn10$data), nrow(svy.asexp.pres08$data), nrow(svy.asexp.pres10$data)))
)
names(svy.asexp.balance.table) <- c("2008", "2010", "2008", "2010")

svy.asexp.balance.caption <- c("This table shows the balance statistics for the Survey 'as an' experiment results. `Previous Turnout' is turnout in the election prior to the survey. `Previous Pres. Party Vote' is vote share of president's party in election prior to the survey. Estimates account for unequal probability of selection and control for precinct size, the variable used to sample precincts into the ``treatment''.")

latex(svy.asexp.balance.table, cdec = c(2,2,2,2), dcolumn = FALSE, ctable = TRUE,
      cgroup = c("Previous Turnout \\%", "Previous Pres. Party Vote \\%"), n.cgroup = c(2,2),
      rowname = c("Estimate", "SE", "95\\% CI", "N"), rowlabel = "",
      col.just = rep("c", 4), insert.bottom = svy.asexp.balance.caption, where = "htp",
      label = "svy.asexp.balance.table", caption = "Balance Statistics Survey for ``as an'' Experiment Results",
      file = "svy.asexp.balance.table.tex")


### Table 8

lm.anychange <- (lm(any.change.pct~tr.list + block, data=exp.data, subset = urban==1))
lm.anychange$se <- vcovHC(lm.anychange,type = "HC0")
lm.anychange.cov <- (lm(paste("any.change.pct~tr.list+block",paste(covariates, collapse="+"), sep = "+"), data=exp.data, subset=urban==1))
lm.anychange.cov$se <- vcovHC(lm.anychange.cov,type = "HC0")

lm.earlyreg <- (lm(earlyreg.changes.pct~tr.list + block, data=exp.data, subset = urban==1))
lm.earlyreg$se <- vcovHC(lm.earlyreg,type = "HC0")
lm.earlyreg.cov <- (lm(paste("earlyreg.changes.pct ~ tr.list + block",paste(covariates, collapse="+"), sep = "+"), data=exp.data, subset=urban==1))
lm.earlyreg.cov$se <- vcovHC(lm.earlyreg.cov,type = "HC0")

lm.latereg <- (lm(latereg.changes.pct~tr.list + block, data=exp.data, subset = urban==1))
lm.latereg$se <- vcovHC(lm.latereg,type = "HC0")
lm.latereg.cov <- (lm(paste("latereg.changes.pct~tr.list+block",paste(covariates, collapse="+"), sep = "+"), data=exp.data, subset=urban==1))
lm.latereg.cov$se <- vcovHC(lm.latereg.cov,type = "HC0")

voterslist.caption <- "The average treatment effect of the complaints treatment on changes to the
voters list. For covariate adjusted results, the adjustment variables
are number of registered voters, 2006 turnout, and 2006 vote share
of the ruling party. Standard errors are heteroskedasticity
consistent and account for blocking."

voterslist_table <- apsrtable(lm.anychange, lm.anychange.cov, lm.earlyreg, lm.earlyreg.cov, lm.latereg, lm.latereg.cov, Sweave = FALSE, omitcoef = expression(grep(pattern=paste("Intercept", paste(covariates,collapse="|"), sep = "|"),x=coefnames)), coef.names = c("Voters List Treatment"), stars = "default", model.names = c("1","2","3","4","5","6"), caption = voterslist.caption, caption.position = "above", label = "tab:dv.voterslist")
voterslist_table[2] <- paste("\n \n  & \\multicolumn{ 2 }{ c }{ Any Changes } & \\multicolumn{ 2 }{ c }{ Early Changes } & \\multicolumn{ 2 }{ c }{ Late Changes }\\\\ \n", voterslist_table[2], sep=" ")
#voterslist_table[1] <- gsub(x = voterslist_table[1], pattern = "\n \n", replacement = "\n \\\\centering \n")
voterslist_table[1] <- gsub(x = voterslist_table[1], pattern = "\\!ht\\]\n", replacement = "\\!ht\\]\n \\\\centering \n")
writeLines(voterslist_table, con = "voterslist_table.tex")


## Table 9

exp.data.labels <- c("Any Change to the Voters List", "Registered Voters (April 17, 2008)", "Pct. Turnout (2006)", "Pct. Vote for President's Party (2006)", "Pct. Registration Changes 1 Week after Intervention", "Pct. Registration Changes 2 Weeks after Intervention", "Pct. Turnout (2008)", "Pct. Vote for President's Party (2008)", "Pct. Reporting Political Pressure by Employers Widespread", "Pct. Reporting Employers Force Employees to Campaign", "Number of Registered Voters (Election Day)")

capture.output(
  stargazer(data.frame(exp.data[, c("any.change.pct", "voters.4.17", "turnout2006", "pctPresParty2006", "earlyreg.changes.pct", "latereg.changes.pct", "turnout08", "pctPresParty08", "pre.job.pressure", "empl.camp", "num.reg.voters"), ]                                                                                                                                                                                                                                                         ),covariate.labels = exp.data.labels, title = "Descriptive Statistics for 2008 Experiment")
  , file = "experiment_descriptive.tex")


## Table 10
capture.output(
stargazer(data.frame(svy10.analysis.data[, c("treatment", "sum.refuse", "sum.dontknow", "UNM.closest", "is.georgia.democracy", "vote.saakshvili.08"), ]), covariate.labels = c("Treatment", "\\# of Refuse to Answer", "\\# of Don't Knows", "Closest to the UMN", "Believes Georgia is a Democracy", "Voted for Saakashvili")), file = "svy_exp_desc.tex")
